# Interpolate points data to raster
library("gstat")
library("ggplot2")
library("tidyverse")
library("automap")
library("fields")

# load data
data <- read.csv("input/data.csv")

# vizualize & evaluate
# factor 1
ggplot(data = data, mapping = aes(x = x, y = y, color = fact1)) + 
  geom_point(size = 2) + 
  scale_color_gradientn(colors = c("blue", "yellow", "red")) + 
  theme_minimal()

# factor 2
ggplot(data = data, mapping = aes(x = x, y = y, color = fact2)) + 
  geom_point(size = 2) + 
  scale_color_gradientn(colors = c("blue", "yellow", "red")) + 
  theme_minimal()


# create grid limits
box <- c(
  "xmin" = min(data$x),
  "xmax" = max(data$x),
  "ymin" = min(data$y),
  "ymax" = max(data$y)
)

# calc step to find 100x100 box grid
boxstep  <- c(
  "x" = as.numeric(box["xmax"] - box["xmin"])/100, ## change 50 to any other grid step, 50 = 50 * 50 cells
  "y" = as.numeric(box["ymax"] - box["ymin"])/100
)

# grid <- expand.grid(x = seq(from = box["xmin"], to = box["xmax"], by = as.numeric(boxstep["x"])),
#                     y = seq(from = box["ymin"], to = box["ymax"], by = as.numeric(boxstep["y"])))

sf_grid <- st_as_sf(data, coords = c("x", "y"))
grid <- sf_grid %>%
  st_bbox() %>%
  st_as_sfc() %>%
  st_make_grid(
    cellsize = c(boxstep["x"], boxstep["y"]),
    what = "centers"
  ) %>%
  st_as_sf() %>%
  cbind(., st_coordinates(.)) %>%
  st_drop_geometry() %>%
  mutate(fact1 = 0, fact2 = 0)

grid_raster <- grid %>%
  rasterFromXYZ()

# plot ALL-IN-ONE picture
par(mfrow=c(2,2))

# fit nearest neighbour interpolation
fit_NN <- gstat(
  formula = fact1 ~ 1,
  data = as(sf_grid, "Spatial"),
  nmax = 10,
  nmin = 3
)

interp_NN <- interpolate(grid_raster, fit_NN)
plot(interp_NN, main = "Nearest neughbour interpolation")

fit_NN.rsm <- sqrt(mean(gstat.cv(fit_NN, debug.level = 0, random = F)$residual^2))

# fit thin plate spline
fit_tps <- fields::Tps(
  x = as.matrix(data[, c("x", "y")]),
  Y = data$fact1,
  miles = F
)

interp_tps <- interpolate(grid_raster, fit_tps)
plot(interp_tps, main = "Thin plate spline interpolation")

# fit kriging ordinary interpolation
fit_krig <- autoKrige(
  formula = fact1 ~ 1,
  input_data = as(sf_grid, "Spatial")
) %>%
  .$krige_output %>%
  as.data.frame() %>%
  dplyr::select(X = x1, Y = x2, Z = var1.pred)

iterp_krig <- rasterFromXYZ(fit_krig)
plot(iterp_krig, main = "Auto Kriging interpolation")

# fit IDW interpolation
fit_idw <- gstat(
  formula = fact1 ~ 1,
  data = as(sf_grid, "Spatial"),
  nmax = 10, nmin = 3,
  set = list(idp = 0.5)
)
interp_idw <- interpolate(grid_raster, fit_idw)
plot(interp_idw, main = "IDW interpolation")

fit_idw.rsm <- sqrt(mean(gstat.cv(fit_idw, debug.level = 0, random = F)$residual^2))
